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ABSTRACT 

In this paper we study an inviscid model for a steady axisymmetric flow 
with swirl. The governing equation is a nonlinear elliptic equation which has 
more than one solution for a certain range of the swirl parameter. The 
physically interesting solutions have closed streamlines that look like vortex 
breakdown ("bubble"-like solutions). A multigrid method is used to find these 
solutions. Using an FMG algorithm (nested iteration), the problem is solved 
in just a few multigrid cycles. 
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1 . INTRODUCTION 


In this paper we study an inviscid model for steady axisymmetric flow with 
swirl, which has solutions with closed streamlines. These solutions have a 
structure similar to that observed experimentally as "bubble"-like solutions 
when vortex breakdown occurs [4]. 

Using a streamfunction-vorticity formulation to the axisymmetric 
incompressible Navier-Stokes equations, it was found [3] that one can reduce 
the problem to a single nonlinear elliptic equation for the streamfunction, in 
case of a special inflow flow and some regularity assumption on the vorticity. 
This nonlinear elliptic equation for the streamfunction has more than one 
solution. The trivial, represents a uniform flow and is of no physical 
interest. The other shows a "bubble"-like structure, the target of our 
numerical study. 

In solving the problem numerically, the problem is reformulated in terras 
of a perturbed streamfunction, i.e., the deviation from the trivial solution. 
In terras of this perturbed streamfunction, the trivial solution is represented 
as an identically zero solution. Our goal then is to find non-zero solutions 
which have "bubble"-like form. 

The approach we have taken in finding these solutions is to seek first for 
a bifurcation point from the trivial branch of solutions. By introducing a 
continuation parameter, we can then start marching on a branch of non-trivial 
solutions that bifurcate from that point. One choice of a continuation 
parameter is arc length [1]. Another choice, which is simpler but may not be 
good in general, is the norm of the perturbed streamfunction. The natural 
parameter in the problem, a swirl velocity parameter, is not good enough since 
it cannot "choose" the non-zero branch as can the former parameter. We 
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therefore choose the norm as a continuation parameter, making the swirl 
velocity parameter an unknown to be determined by the solution. 

The multigrid approach used for solving the problem is similar to the one 
used in [5] for solving the Bratu problem. The relaxation in this method 
consists of three steps: (i) a local relaxation to smooth the error; (ii) a 
step to update the norm of the solution; and (iii) a step to update the swirl 
velocity parameter. An FMG algorithm (nested iteration) is used. That is, a 
solution for the prescribed norm is found first on the coarsest level, and 
then interpolated to finer levels, where on each level a few basic multigrid 
V-cycles are performed before proceeding to yet finer level. 

The coarsest level, when solved to get an initial approximation for finer 
levels, uses a continuation method. Here the problem was solved first for a 
small norm, and then the norm is gradually increased until the prescribed norm 
is reached. Each time the norm is increased, the solution of the previous 
step was used as initial approximation. By solving for a bifurcation point 
from the trivial solution, a first approximation for the smallest norm problem 
was obtained. 

Once a solution on the coarsest level is obtained for a prescribed norm, 
it is possible to solve finer grid problems without continuation. 

The same problem we are discussing here was treated by a completely 
different method and is reported in [3], There, a single grid method was used 
with a least squares formulation of the problem. The amount of work needed 
for that approach is considerably larger than the one reported here. Computed 
solutions by the two different formulations are in good agreement. 
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2 . OH DERIVATION OF THE GOVERNING EQUATION 

We summarize here the derivation of the equations used in the numerical 
process as given in ( 3 ], In cylindrical coordinates (x,r, 0 ) the 
incompressible Navier-Stokes equations can be written in terms of a stream- 
function tji, vorticity m, and circulation k as 


(ml + Cm.) + 3 = Ite [“rr + 7 “r 

r x L 

uk + wk = 4r* |"k - — k +k 1 

r x Re |_ rr r r xxj 




(2.1b) 

(2.1c) 


where k = rv, w = w^ — u^ and Re is the Reynolds number. The velocity 
components in the x,r,0 directions are w, u, v, respectively, of which w 
and u are given in terms of the streamfunction by 


w 



(2.2a) 


u = 



(2.2b) 


It is shown in [3] that in the inviscid case (Re = ~), one finds that the 
circulation k and the vorticity u are functions of the streamfunction i|» 
only. Therefore, k and a) can be determined outside the "bubble" from the 
inflow boundary condition. In the model discussed it is assumed that the same 
functional dependence of k, w on is true also inside the bubble 
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(negative ty). This imposes some regularity on the solution* 
For the inflow conditions 


|V 0 r(2 - r 2 ) r < 1 

v(0,r) = ] 

( V 0 /r r > 1, 


(2.3a) 


w(0,r) = 1, 


(2.3b) 


it is possible to write k and u in terms of the streamfunction 


as 


k (0, r) = 


16 V 2 ^ 2 (1 


- 


w(0,r) = 


16 V 2 (l + 2^ 2 - 3tJ/)(r 2 /2 - tp) 


0 


^ < V2 
> V2 

'l' < V2 
'P > V2 


(2.4a) 


(2.4b) 


and therefore, the equation obtained for is 


r <V r) r + ^ xx " “ 4v q " r Z /2) 


(2.5a) 


where 


a 2 (<P) = 


’4(1 + 2tp - 3tp) 
0 


tp < V2 

> V2 


(25b) 


The reduction of the governing equations to a single nonlinear elliptic 
equation is possible if the relation tp = f(r) in the inflow boundary can be 
inverted to get r = g(tp). When g( xp) is introduced in the expression for 
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v at the inflow boundary, one has v as a function of \|) in that boundary 

and therefore k(t}>)> Note that, in general, one cannot expect to 

analytically invert the relation = f(r), and so the reduction of the 

governing equations is possible only for very special inflows. 

2 

Numerical experiments were done in terms of <J> = <J> - ■=— , which is a 

2 

perturbation from the trivial solution 4* = — that represents a uniform flow. 


3. NUMERICAL ALGORITHM 
3.1. Discretization 

2 

The equation for <f> = 4* - r /2 is given by 

rf- ♦ ) + ♦ + 4 V 2 a 2 (<f>)<p = 0, R » (0,a) x (0,b) (3.1a) 

v r T r'r xx 0 

<*> = 0 , on 3 ft (3.1b) 


where 


a 2 (*) 


4 $ - 1 + I- (2+ - 1 + r ) 


otherwise 


Equations (3.1) are discretized as 


’ 1+1,3 24> ij + *1-1 ,j fj [ 2 , h 

h 2 h 2 Vi + r j iJ+1 


t 2 2 f , h \ ,h 

f Q “ (<i> ij^ij = °» in n 


(3.1c) 


(3.2a) 


- 0, on 3R h 


(3.2b) 
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where fl h = {nh.rah), 0 < nh < a, 0 < mh < b}. 


3 »2 . General Strategy for Solving the Discretized Equations 

Equation (3*2) has the trivial solution <f>^ = 0 for any Vq. This 

solution corresponds to a uniform flow and is not interesting physically. We 
seek solutions which represent vortex breakdown so that II <f> ^ II ^ * 0, where 

h 2 2 2 

II4-V = h z Z (3.3) 

Iterating on equation (3.2) by any iterative method may lead us to the trivial 
solution. In order to rule out this possibility, we specify the norm of the 
discrete solution we want to find, while making free the swirl velocity 
parameter Vq. 

To summarize, we solve equation (3.2) for (<}> h , V Q ) under the 
constraint 

"<f> h|12 " g 0 , (3.4) 

where gQ is given. 

A relaxation scheme for (<f> h ,V 0 ) in equation (3.2) together with the 
constraint (3.4) is described next. 


3.3. Relaxation 

Equations (3.2), (3.4) form a nonlinear system of equations for (^Vq). 
The relaxation used for this system has three steps: (i) a local process for 
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sraoothing <}>^ in equation (3.2); (ii) a global change to satisfy (3.4); and 
(iii) updating the swirl parameter Vq. That is, one relaxation consists of 
doing (i), (ii), and (iii) successively. 

(i) local relaxation 

Scan the point (i,j) € fl in lexicographic ordering; at each point 
(i,j) solve (3.2) approximately for by applying one Newton 

iteration. 

(ii) global step 

Compute 3 = / gp/ll^ll 2 . 

Then make the change 


♦Jj * 6 +« ’ 


(i,j) € fl. 


(iii) updating Vq 

Change Vq such that the following equation holds 


<L h <{> h + 4 Vq a 2 (<}> h )<j> h , <j> h > = <f h ,<J> h > 


(3.5) 


H Vi 1 

where L $ is the discretization of L<f> = r(-— <p ) + <J> , <•,•> 

r r r xx 

2 i 

denotes the inner product, <u,v> = h J u.. v.., and f n is the 

ij 1J 13 

right-hand side of equation (3.2). (In a multigrid process f is 
nonzero on coarse grids.) 


We now come to the description of the multigrid algorithm used to solve (3.2), 
(3.4) for (<j, h , V Q ). 
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3.4.1. Basic Cycle : 

Given a sequence of discretizations with mesh sizes 

h r > h 2 >...> h m , where h k = 2h k+1 . The h k -grid equation is generally 
written as 

T k .k c k 

L * = f (3.6) 

where L k approximates L k+1 (k < m) (e.g., they all are finite-difference 
approximations to the same differential operator). The algorithm for 
improving a given approximate solution $ k to (3.6) is denoted by 

? k * MG(k, $ k , f k ) (3.7) 

and is defined recursively as follows: 

If k = 1, solve (3.6) by several relaxation sweeps; qtherwise do steps 

(A) - (D): 

(A) Perform Vj relaxation sweeps on (3.6), resulting in a new 
approximation <{> . 

(B) Starting with ^ = I k ^ '$ k , perform one cycle 

? k - 1 . t- 1 , L k_1 ? k -‘ + ^-‘(f k - L k **)). 

(C) Calculate 

* k = t + i k _i^ k_1 - lk_1 

(D) Perform v 2 additional relaxation sweeps on (3.6) starting with 

— i 

<J> and yielding the final <j> of (3.7). 
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In this algorithm I k_1 , T^ -1 are f ine-to-coarse grid transfer operators; 
I k 1 is an interpolation operator. We refer to the above cycle as MG(v 1> v 2 ). 
In the notation of this section (3.6) includes both equations (3.2) and (3. A). 

The basic cycle described above is for improving a given approximation on 
level k. The full multigrid (FMG) process involves solving the problem on 
the coarsest grid, interpolating it to finer grids, and making the cycle 
MG(Vj,v 2 ) a few times after each refinement. 


3. A. 2. Full Multigrid Algorithm (FMG) 

1. Solve (3.6) for k = 1, using a continuation method (see remark 
below). 

2. Set k = k + 1 and 

= n k _^ ^ , where JI k _^ is a bicubic interpolation. 

3. Perform y(k) times the cycle 
? k «■ MG(k, * k , f k ). 

A. If k < m, go to step 2; otherwise stop. 

A Remark on Step 1 of the FMG Algorithm (Continuation Method) 

Since the problem involved is a nonlinear one, and we are using a Newton 

iteration, a good initial approximation may be needed to get fast convergence 

for k = 1 (the coarsest grid). This has been achieved by using a continuation 

h 2 

process where we solve first for a small norm 11$ H > then gradually 
increasing it until the prescribed norm is obtained. Each time the norm is 
increased, the solution of the previous step is used as an initial 
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approximation. In order to get a good initial approximation for the smallest- 
norra problem, we have solved for the bifurcation point from the trivial branch 
of solutions. 


3 *5 Solving for the Bifurcation Point 

At a bifurcation point (<J> , V^), the linearized problem of (3.1) must 
have a zero eigenvalue, and the corresponding eigenfunction gives rise to a 
second branch of solutions. Since <{> = 0 is a solution for any Vq, we may 
try to find a bifurcating branch from the trivial one (0,V Q ). The 
linearized equations around (0 ,Vq) are given by 

W xx + r ^7 W Jr + 4V 0 “ 2 <°) w = °» in ft (3.8a) 

W = 0, on an. (3.8b) 

If there exists a bifurcating branch from the trivial one (0,V Q ), equation 

(3.8) has a solution (W ,V Q ) with HW » 2 = 1 where n n 2 denotes the L 2 
norm. 

We discretize (3.8) in a way similar to the discretization of (3.1). The 
constraint 

nw h n 2 = i. 


is added to ensure a non— zero solution to the problem. The process of solving 
the eigenvalue problem is identical to the process of solving (3.2), (3.4). 



- 11 - 


Once this linear eigenvalue problem is solved, we can use = ±eW as 
an initial approximation for our original problem with a prescribed norm of 
e. The sign is chosen such that <|>q has negative values, to ensure that the 
total streamfunction ~ “ + ♦ will have closed streamlines with negative 
values (the bubble). 


4. NUMERICAL RESULTS 

Experiments were performed with equations (3.2), (3.3) using FMG 
algorithm of Section 3.4.2. In these experiments the domain was 

J2 h = {(nh, £h) , 0 < nh < 5, 0 < £h < 2} . 

Three levels were used in the multigrid algorithm where the finest grid 

problem has mesh size 1/16. On the coarsest level 20 relaxations were 

performed while on finer grids = 3, y(^) =4. In all numerical 

k-1 — k-1 k-1 

experiments 1^ = 1^ is injection, 1^ is bilinear interpolation, and 

n k-1 is bicubic interpolation. 

2 

Tables I-IX contain the l^-norra of the residuals and the values of Vq 
at the end of each cycle on the finest grid. Cycle #0 refers to the 
approximation obtained from the previous level as an initial guess. Figures 
1-9 show the streamlines (contours of tj») for the different cases. The value 

A tfc 

of Vq , the swirl parameter value for which bifurcation occurs is Vq = 1.0069 
(computed on coarsest level). 

The experiments clearly show that the multigrid method suggested is very 

2 

efficient. In fact, as seen by the convergence history for Vq, it is enough 
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to take y(k) - 2, Instead of y(k) = 4, i.e., by 2 FMG cycles the problem is 
already solved. 

The results show that bigger bubbles are obtained for smaller swirl 
parameters, contradicting to what one would expect. This may be the result of 
the assumption made in the model, that the same functional dependence of 
k, u on i|i holds Inside as well as outside the bubble. A future study will 
investigate this point by solving the full systems (2.1), making no extra 
assumptions. 
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Table I. H(J> h II 2 = .005 


cycle it 

II Residuals B 2 

CM O 
> 

0 

.362 (-1) 

.95088 

1 

.986 (-3) 

.96069 

2 

.843 (-4) 

.96039 

3 

.148 (-4) 

.96041 

4 

.745 (-5) 

.96042 


Table 

III. 

it4> h 

2 

ir = ,ii 


cycle 

it 

II Residuals II £ 

CM O 
> 

0 


.122 


.59214 

1 


.233 

(-2) 

.54739 

2 


.215 

(-3) 

.54732 

3 


.615 

(-4) 

.54733 

4 


.542 

(-4) 

.54733 


Table V. 

ii 4> h n 2 = .2 


cycle if 

II Residuals H 2 

< 

0 

.150 

.42902 

1 

.242 (-2) 

.43301 

2 

.193 (-3) 

.43294 

3 

.366 (-4) 

.43294 

4 

.266 (-4) 

.43294 


Table 

II. 

04 > h n 2 = .05 


cycle 

if 

HResiduals If 2 

V 2 

0 

0 


.948 (-1) 

.68322 

1 


.232 (-2) 

.68962 

2 


.251 (-3) 

.68939 

3 


.113 (-3) 

.68941 

4 


.918 (-4) 

.68941 

Table 

IV. 

n<f> h n 2 = .15 


cycle 

if 

II Residuals II 2 

V 2 

v o 

0 


.135 

.48347 

1 


.243 (-2) 

.48803 

2 


.168 (-3) 

.48798 

3 


.474 (-4) 

.48798 

4 


.425 (-4) 

.48798 

Table VI. 

ll<J> h ll 2 = .4 


cycle 

if 

II Residuals II 2 

'2 

0 


.192 

.30435 

1 


.271 (-2) 

.30725 

2 


.239 (-3) 

.30719 

3 


.177 (-3) 

.30719 

4 


.176 (-3) 

.30719 
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Table VII. n<J> h n 2 = .6 


cycle # 

0 Residuals II 2 

CM O 
> 

0 

.230 

.24006 

1 

.303 (-2) 

.24335 

2 

.218 (-3) 

.24231 

3 

.188 (-3) 

.24231 

4 

.175 (-3) 

.24231 


Table 

IX. 

Il4> h n 2 

= 2.0 


cycle 

# 

11 Residuals H 2 

CM O 
> 

0 


.428 


.10176 

1 


.701 

(-2) 

.10276 

2 


.777 

(-3) 

.10275 

3 


.584 

(-3) 

.10275 

4 


.574 

(-3) 

.10275 


Table VIII. tl<}> h ll 2 = 1.0 


cycle # 

II Residuals II 2 

CM O 
> 

0 

.295 

.17139 

1 

.385 (-2) 

.17303 

2 

.363 (-3) 

.17302 

3 

.294 (-3) 

.17302 

4 

.278 (-3) 

.17302 
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.48798. 
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